% 本例程用于通过子载波系数增加线性相位的方式
% 产生过采样的时域序列

clc;close all;
% 子载波系数数量
N=128;
% 过采样率
os=4;

% 产生子载波系数，该系数在无过采样时为恒模序列
% x=fft(exp(1i*2*pi*rand(N,1)));
x = exp(1i*2*pi*rand(N,1));

% 第一组采样序列
y1 = ifft(x);
% 第二组采样序列
x1 = x .* exp(1i*2*pi/(os*N)*(0:N-1)');
y2 = ifft(x1);
% 第三组采样序列
x1 = x1 .* exp(1i*2*pi/(os*N)*(0:N-1)');
y3 = ifft(x1);
% 第四组采样序列
x1 = x1 .* exp(1i*2*pi/(os*N)*(0:N-1)');
y4 = ifft(x1);

% 最终采样序列
y = [y1,y2,y3,y4];
y = y.';
y = y(:);

% 通过ifft产生过采样序列
yos = ifft(x, N*os);

figure(1);hold on;
subplot(211);hold on;title("频谱+线性相位的方式获得过采样序列");
stem((1:os:os*N),abs(y1));
stem((2:os:os*N),abs(y2));
stem((3:os:os*N),abs(y3));
stem((4:os:os*N),abs(y4));
subplot(212);hold on;title("直接IDFT产生过采样序列")
plot(abs(yos))

% figure(2);hold on;
% linePha = exp(1i*2*pi*1/(os*N)*(0:N-1)');
% temLinePha = ifft(linePha);
% subplot(211);hold on;title("线性相位时域幅度特性")
% stem(abs(temLinePha))
% subplot(212);hold on;title("线性相位时域相位特性")
% plot(angle(temLinePha))

% 以下例程研究对比了不同恒包络序列DFT后自相关的情况

load('smoth_weights_pca_l1.mat', 'smoth_weights');
N = length(smoth_weights);
smothSeq = ifft(smoth_weights, os*N);
smothSeq = sqrt(os*N) * smothSeq / norm(smothSeq, 2);
fftOfUmseq = fft(smothSeq);
% fftOfUmseq = fft(exp(1i*2*pi*rand(os*N,1)));
convOfUmseqDFT = periodic_autocorr(fftOfUmseq);
figure(3);hold on;title("时域恒包络与频谱自相关的关系")
legend();
plot(abs(convOfUmseqDFT), 'DisplayName', '平滑系数周期自相关');
plot(os*N*abs(fft(ones(os*N,1))), '--', 'DisplayName', '全1序列DFT')
temp = fft(ones(os*N,1));
xortemp = periodic_autocorr(temp);
plot(abs(xortemp), '-.', 'DisplayName', '全1序列系数自相关');

load('non_smoth_weights.mat', 'non_smoth_weights');
N = length(non_smoth_weights);
nonSmothSeq = ifft(non_smoth_weights, os*N);
nonSmothSeq = sqrt(os*N) * nonSmothSeq / norm(nonSmothSeq, 2);
fftOfUmseq = fft(nonSmothSeq);
% fftOfUmseq = fft(exp(1i*2*pi*rand(os*N,1)));
convOfUmseqDFT = periodic_autocorr(fftOfUmseq);
plot(abs(convOfUmseqDFT), 'DisplayName', '非平滑系数周期自相关');

function R = periodic_autocorr(x)
    N = length(x);
    R = zeros(1, N);
    for k = 0:N-1
        for n = 0:N-1
            idx = mod(n + k, N) + 1;  % MATLAB 下标从 1 开始
            R(k + 1) = R(k + 1) + x(n + 1) * conj(x(idx));
        end
    end
end


